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1 Introduction 

Efforts to develop better and more efficient algorithms for Monte Carlo simula- 
tions have a long history, in which the Fortuin-Kasteleyn (FK) transformation 
[0 for Potts models has played a pivotal role. The most common use of this 
transformation has been to create algorithms in which clusters of spins in the 
Potts model are flipped simultaneously j^, [|. In this paper, we present a new 
algorithm using the FK representation in the spirit of recent work by Newman 
and Ziff on the percolation problem. The method produces independent 
samples and sums up a large number of configurations for each sweep. The par- 
tition function and thermodynamic averages for all values of the temperature T 
can be computed from a single run. 

Consider the Potts model with the Hamiltonian, 



where the summation is over nearest neighbor pair, Cj = 1,2,..., q. The FK 
transformation allows us to write the partition function in the percolation rep- 
resentation as 



where we sum over all configurations of bonds connecting nearest neighbor sites, 
p = 1 — exp(— J/(kT)Y b is the number of bonds present, M — dN is the 
maximum number of possible bonds, and N C (T) is the total number of clusters 
for a given configuration of bonds. 

In order to evaluate thermodynamic averages in this representation, we carry 
out the summation over all configurations in two steps. First, we sum over the 
number of bonds 6, and then for each value of 6, we sum over all configurations 
consistent with that number of bonds. Thus we can write 
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where Tb is a configuration with exactly b bonds. If wc could compute Cb and 
Qb for some observable Q(T) for every b, then we could compute the function 
at any values of p or T. 

The problem then reduces to the computation of the normalization constants 
Cb and the expectation values in an ensemble with probability distribution pro- 
portional to q N °( r \ For the special case of q = 1 (bond percolation), the values 
of Cb are the binomial coefficients and the configurations with fixed number 
of bonds are weighted uniformly. For general case, we describe two sampling 
methods. The first is very simple, and useful conceptually, but exponentially 
inefficient. The second turns out to be very efficient. 



2 A Survival and Death Process 

Starting with an empty lattice, one sweep consists of repeated application of 
the following steps until the process dies: 

1. Pick an unoccupied neighbor pair at random for the next bond. 

2. If inserting a bond 

(a) does not change the cluster number (A7V C = 0), accept the configu- 
ration; 

(b) merge two clusters, so that the cluster number decreases by 1, (AN C = 
— 1), accept the configuration with probability 1/q, or reject the con- 
figuration and terminate the process (and begin the next sweep from 
an empty lattice). 

3. Take statistics of the survival configurations (with equal weights). 

The probability distribution of the configurations generated in this manner 
with b bonds is 

P(T b )= b -^^q^-^. (5) 

The survival probability after b bonds is then 
( . b\(M-b)\ _ N 

r b 

where N is the total number of sites. As the survival probability decays expo- 
nentially with the number of bonds added, we expect very poor statistics for 
large values of b. However, the eventual survival probability, Sm, is equal to 
q~ N+1 , which is nonzero. 
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3 Binary Tree Summation Method 



To deal with the exponentially decreasing survival probability, we have devel- 
oped the following method for summing multiple bond sequences in a single 
sweep. For a configuration Tf, we count the number of type-0 bonds hq of un- 
occupied pairs of sites that are on the same cluster, and the number of type-1 
bonds n\ of unoccupied pairs of sites that would connect two different clusters. 
Clearly no + ii± + b = M. At each step, we pick a type-1 bond with equal prob- 
ability from among all current type-1 bonds. Starting with an empty lattice 
and merging clusters at each step, we continue until all sites are members of 
the same cluster. Along the way, we collect statistics Q(i) for each of the N 
configurations generated. 

After a sweep has been carried out, we construct all possible paths that 
a full simulation of both type-0 and type-1 bonds would have taken if we had 
joined clusters in exactly the same sequence as in the actual simulation, but had 
also inserted a random number of type-0 bonds at each step. At each choice 
between types of bond, the type-1 bond is given a relative weight of n\j q and 
a type-0 bond rig. The total weight for a path starting from an empty lattice 
to a particular configuration is the product of the factors uq or ri\jq depending 
on the path taken. 

To understand this algorithm, imagine that we had actually followed the 
branching process. At each step, a configuration may split into two configura- 
tions, one with a type-1 bond added (with probability 1/rii) and one with a 
type-0 bond added (with probability of 1/no). The probability of appearance 
of a particular path is ri[ n /(fc)(rfc)] ; where f(k) = or 1 depending on the 
choice of type-0 or 1 bond. In taking statistics, we have weighted with the 
inverse of the factor, multiplied by additional factors of q proportional to the 
number of clusters N C (T). The net effect is the required sample average with 
overall weight g^T) _ Now the key observation is that we do not need an explicit 
simulation for the type-0 bonds, since the type-0 bonds have no effect on the 
measured quantities. 

The binary tree summation algorithm has several attractive features. First, 
the usual slowing down due to correlation between samples is absent. Each 
sweep is independent. Second, for each sweep, data for a very large number of 
samples are collected. Although they are highly correlated, an exponentially 
large number of paths can be summed efficiently with 0(N 2 ) operations per 
sweep. Third, unlike multicanonical simulations || or the flat histogram method 
[|J, there are no unknown weighting factors to determine. 

4 Implementation 

Let w(b, i) be the weight of the total contributions from all possible paths to 
the state specified by the number of bonds b and merge sequence number i. 
These quantities can be calculated from the starting condition w(0, i) = t^o, 
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the constraint that no(b,i) > 0, and the recursive equation 

w(b + 1, i) = w(b, i)n (b, i) + w(b, i — l)rii(&, i — l)/q, (7) 

where n\(b,i) = n\{i), no(b,i) — no(i) — b + i. The value w(b,i) is nonzero only 
for b > i. The computation of the weights w(b, i) is similar to that of binomial 
coefficients. The final contribution to the statistics at b number of bonds is 

JV-l 

Q b =(W b )- 1 (Y,™(b,i)Q(i)), (8) 

i=0 

where the average is over simulation sweeps, and W b = w(b, i) is the total 
weight at a given b. We can then show that 

<"« " W^ s >- (9) 

which allows us to compute c&. The result is numerically identical to computing 
the conditional survival probability Sb+i/Sb from the expectation value of Q = 
(n + n 1 /q)/(M-b). 

During the simulation, the values of no and n\ can be updated efficiently. 
For each cluster, we keep a list of unoccupied bonds with other clusters. When 
two clusters (A and B) are joined, we merge the smaller one with the larger 
one, and remove and count the number of bonds tiab connecting the clusters. 
We update according to no <— uq + tiab — 1, rii <— n-i — tiab- The timing of 
our program shows that this part of the algorithm scales nearly linearly with 
number of sites TV, as expected. 



5 Results and Discussions 



We note that the coefficients q, are related to the density of states n(E), 
which gives the coefficients of the partition function polynomial in the variable 
expf— J/(kT)J. In the FK percolation representation, it becomes a polynomial 
in —p). By a proper change of variables, we can find exact results of q, 
for the two-dimensional Ising model. 

We define the following errors to test our method against exact results for 
the two-dimensional Ising model: 
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Due to a special cancellation for this algorithm, eo is exactly zero. The error ei, 
along with the maximum and average errors in the energy and specific heat, are 
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Table 1: Errors for 10 6 sweeps with respect to the exact results of an L x L Ising 
model. While t\ is taken from an average over many runs, the other results are 
from a single run. The CPU times (on a 1.53GHz Athlon) are in units of 10 -6 
second per sweep per site. 



LA 8 16 32 50 



cpu t 


1.88 


2.81 


6.92 


24.0 


72.6 


fl 


0.0000634 


0.000178 


0.00049 


0.0032 


0.031 


'max 


0.000128 


0.000113 


0.00012 


0.0015 


0.0068 


fc AVE 


0.0000184 


0.0000103 


0.0000055 


0.000046 


0.00014 


e MAX 


0.000306 


0.00046 


0.00124 


0.0177 


0.096 


e AVE 


0.000034 


0.000031 


0.000040 


0.00050 


0.0020 



listed in Table 1. Since the algorithm asymptotically takes 0(N 2 ) operations 
per sweep, a fair comparison with other methods should compare the total CPU 
times. A comparable N-fold way transition matrix Monte Carlo (TMMC) run 
took 1.9 microsecond per sweep per site on the same machine. Thus, the present 
method is superior for small lattices with linear size L < 16. For L = 32 it is 
comparable to TMMC [[|. For much larger lattices, it becomes less favorable 
mainly due to the 0(N 2 ) nature of the algorithm. It is quite likely that we can 
speed up the computation using special properties of the weights. 

Our method is applicable for Potts models with any number of states, in- 
cluding fractional or negative values. Work is currently in progress to apply this 
algorithm to a number of problems of interest. 
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